Environmental effects on corticosterone in nestling tree swallows

Last updated: November 17, 2025

Abstract:

For an individual, the ability to cope with challenges is fundamental to surviving in a changing environment. The corticosterone stress response has a major impact on the ability of an organism to react to the challenges of rapid environmental change. Previous research has shown that when faced with a prolonged cold snap, adult female tree swallows increase their baseline and stress-induced corticosterone levels, and that individual variation in corticosterone levels predicts reproductive success during challenges. Tree swallows are advancing reproduction with warming springs and are more likely to encounter cold snaps during breeding, which can reduce reproductive success. Here, using data from a long-term study of free-living tree swallows collected from 2013-2025, we show that nestling tree swallows show an immediate increase in corticosterone levels when ambient temperatures fall below a threshold of 18.5°C, which has previously been identified as the threshold at which insect availability begins to decline. As a result of earlier reproduction, nestlings are more likely to experience cold snaps during development, making it especially important to understand their effects. Ongoing work is also testing whether male and female nestlings display dimorphism in their corticosterone responses to further our understanding of how nestlings are responding to thermal environmental shifts.

Research Questions:

How does temperature influence nestling corticosterone levels?
Do male and female nestlings differ in their corticosterone regulation?
Does sex mediate the effect of temperature on the corticosterone response in nestlings?

Hypotheses:

Lower temperatures during developmental periods lead to increased corticosterone responses;
Corticosterone regulation in nestlings will be mediated by sex; and
lower temperatures during developmental periods induce sex-specific effects on the stress response in nestling tree swallows.

Predictions:

Hypothesis 1: Temperature will show a negative relationship with baseline and stress-induced corticosterone, with there being an increase in these measures when both average daily temperatures and average 3-day temperatures are lower; Temperature will have no relationship with post-dex corticosterone levels.

Hypothesis 2: Males and females will show base differences in their corticosterone regulation, with females showing lower corticosterone across all measures, based on previous studies done in adult birds (baseline corticosterone, stress-induced corticosterone, and post-dex corticosterone).

Hypothesis 3:Females will showcase a stronger corticosterone response (across all measures) in response to lower temperatures, consistent with prior work in our population showing that females may more sensitive to some stressors.

Question 1: How does temperature influence nestling corticosterone levels?

In this section, I am running GAMM models, as based on prior checking, the effect of temperature on corticosterone is non-linear, and a smooth term fits the data better than a linear term. Average daily temperature is used as the predictor, as this seemed to have a slightly stronger effect than average 3 day temperature, and all corticosterone measures have been log-transformed to conform to a normal distribution better.

Baseline corticosterone

bleed1_weather_gamm4 <- gamm4( # run gamm model 
  log_bleed1_final_cort ~ s(avgC_capture_day) + (bleed1_latency_sec),
  random = ~(1 | nest_key),
  data = weather_df
)
summary(bleed1_weather_gamm4$gam) # view results of gam

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log_bleed1_final_cort ~ s(avgC_capture_day) + (bleed1_latency_sec)
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.993948   0.114410   8.688   <2e-16 ***
## bleed1_latency_sec 0.001967   0.001082   1.818   0.0695 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df     F p-value    
## s(avgC_capture_day) 3.703  3.703 23.39  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.206   
## lmer.REML =   1709  Scale est. = 0.73143   n = 601

summary(bleed1_weather_gamm4$mer) # view random effects

## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 1709
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8056 -0.6282 -0.0074  0.5488  3.5072 
## 
## Random effects:
##  Groups   Name                Variance Std.Dev.
##  nest_key (Intercept)         0.3886   0.6233  
##  Xr       s(avgC_capture_day) 1.2794   1.1311  
##  Residual                     0.7314   0.8552  
## Number of obs: 601, groups:  nest_key, 164; Xr, 8
## 
## Fixed effects:
##                          Estimate Std. Error t value
## X(Intercept)             0.993948   0.113681   8.743
## Xbleed1_latency_sec      0.001967   0.001082   1.818
## Xs(avgC_capture_day)Fx1 -0.556697   0.323104  -1.723
## 
## Correlation of Fixed Effects:
##             X(Int) Xbl1__
## Xbld1_ltnc_ -0.841       
## Xs(vgC__)F1 -0.029  0.032

plot(bleed1_weather_gamm4$gam, pages = 1) # view smooth term

k.check(bleed1_weather_gamm4$gam) # check the k to see if the model is "wiggly" enough

##                     k'      edf   k-index p-value
## s(avgC_capture_day)  9 3.703261 0.8236044       0

par(mfrow = c(2, 2))
gam.check(bleed1_weather_gamm4$gam) # check residuals to make sure model is normal

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                      k' edf k-index p-value    
## s(avgC_capture_day) 9.0 3.7    0.82  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

base_cort <- draw(bleed1_weather_gamm4$gam, term = "avgC_capture_day")$data %>%
  dplyr::select(avgC_capture_day = avgC_capture_day, fit = .estimate, se = .se, cil = .lower_ci, ciu = .upper_ci)

base_cort_plot <- ggplot() +
  geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
  geom_line(data = base_cort, aes(x = avgC_capture_day, y = fit), color = base_color, linetype = "solid", size = 0.3) +
  geom_ribbon(data = base_cort, aes(x = avgC_capture_day, ymin = cil, ymax = ciu), fill = base_color, alpha = 0.4) +
   labs(x = NULL, y = "Predicted Partial Baseline Corticosterone") +
  theme_classic() +
  scale_x_continuous(breaks = seq(0, 30, 10)) +
  xlim(c(13, 30)) +
  theme(panel.grid = element_blank(), axis.text = element_text(size = 12),axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

base_cort_plot

Stress-induced corticosterone

bleed2_weather_gamm4 <- gamm4( # run gamm model 
  log_bleed2_final_cort ~ s(avgC_capture_day) + (log_bleed1_final_cort),
  random = ~(1 | nest_key),
  data = weather_df
)
summary(bleed2_weather_gamm4$gam) # view results of gam

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log_bleed2_final_cort ~ s(avgC_capture_day) + (log_bleed1_final_cort)
## 
## Parametric coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.32155    0.07382  31.447  < 2e-16 ***
## log_bleed1_final_cort  0.22625    0.03584   6.313 5.54e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df     F  p-value    
## s(avgC_capture_day) 3.544  3.544 6.734 0.000112 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.219   
## lmer.REML = 1447.7  Scale est. = 0.48454   n = 572

summary(bleed2_weather_gamm4$mer) # view random effects

## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 1447.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1550 -0.5149  0.0832  0.5904  2.3502 
## 
## Random effects:
##  Groups   Name                Variance Std.Dev.
##  nest_key (Intercept)         0.4551   0.6746  
##  Xr       s(avgC_capture_day) 1.2671   1.1256  
##  Residual                     0.4845   0.6961  
## Number of obs: 572, groups:  nest_key, 158; Xr, 8
## 
## Fixed effects:
##                         Estimate Std. Error t value
## X(Intercept)             2.32155    0.07431  31.242
## Xlog_bleed1_final_cort   0.22625    0.03606   6.275
## Xs(avgC_capture_day)Fx1 -0.45712    0.32667  -1.399
## 
## Correlation of Fixed Effects:
##             X(Int) Xl_1__
## Xlg_bld1_f_ -0.542       
## Xs(vgC__)F1 -0.038  0.069

plot(bleed2_weather_gamm4$gam, pages = 1) # view smooth term

k.check(bleed2_weather_gamm4$gam) # check the k to see if the model is "wiggly" enough

##                     k'      edf   k-index p-value
## s(avgC_capture_day)  9 3.544486 0.6504017       0

par(mfrow = c(2, 2))
gam.check(bleed2_weather_gamm4$gam) # check residuals to make sure model is normal

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                       k'  edf k-index p-value    
## s(avgC_capture_day) 9.00 3.54    0.65  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

stress_cort <- draw(bleed2_weather_gamm4$gam, term = "avgC_capture_day")$data %>%
  dplyr::select(avgC_capture_day = avgC_capture_day, fit = .estimate, se = .se, cil = .lower_ci, ciu = .upper_ci)

stress_cort_plot <- ggplot() +
  geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
  geom_line(data = stress_cort, aes(x = avgC_capture_day, y = fit), color = purple_color, linetype = "solid", size = 0.3) +
  geom_ribbon(data = stress_cort, aes(x = avgC_capture_day, ymin = cil, ymax = ciu), fill = purple_color, alpha = 0.4) +
  labs(x = NULL, y = "Predicted Partial Stress-Induced Corticosterone") +
  theme_classic() +
  scale_x_continuous(breaks = seq(0, 30, 10)) +
  xlim(c(13, 30)) +
  theme(panel.grid = element_blank(), axis.text = element_text(size = 12), axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15))

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

stress_cort_plot

Dex corticosterone

bleed3_weather_gamm4 <- gamm4( # run gamm model 
  log_bleed3_final_cort ~ s(avgC_capture_day) + (log_bleed1_final_cort),
  random = ~(1 | nest_key),
  data = weather_df
)
summary(bleed3_weather_gamm4$gam) # view results of gam

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log_bleed3_final_cort ~ s(avgC_capture_day) + (log_bleed1_final_cort)
## 
## Parametric coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.77662    0.06524  27.233   <2e-16 ***
## log_bleed1_final_cort -0.01010    0.03350  -0.301    0.763    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df     F p-value   
## s(avgC_capture_day) 2.852  2.852 4.569 0.00815 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0417   
## lmer.REML = 954.39  Scale est. = 0.32079   n = 454

summary(bleed3_weather_gamm4$mer) # view random effects

## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 954.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5213 -0.5110 -0.0397  0.5311  5.1130 
## 
## Random effects:
##  Groups   Name                Variance Std.Dev.
##  nest_key (Intercept)         0.2768   0.5261  
##  Xr       s(avgC_capture_day) 0.3262   0.5711  
##  Residual                     0.3208   0.5664  
## Number of obs: 454, groups:  nest_key, 121; Xr, 8
## 
## Fixed effects:
##                         Estimate Std. Error t value
## X(Intercept)             1.77662    0.06500  27.333
## Xlog_bleed1_final_cort  -0.01010    0.03337  -0.303
## Xs(avgC_capture_day)Fx1  0.24314    0.20741   1.172
## 
## Correlation of Fixed Effects:
##             X(Int) Xl_1__
## Xlg_bld1_f_ -0.516       
## Xs(vgC__)F1  0.047 -0.096

plot(bleed3_weather_gamm4$gam, pages = 1) # view smooth term

k.check(bleed3_weather_gamm4$gam) # check the k to see if the model is "wiggly" enough

##                     k'      edf   k-index p-value
## s(avgC_capture_day)  9 2.851564 0.7074518       0

par(mfrow = c(2, 2))
gam.check(bleed3_weather_gamm4$gam) # check residuals to make sure model is normal

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                       k'  edf k-index p-value    
## s(avgC_capture_day) 9.00 2.85    0.71  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dex_cort <- draw(bleed3_weather_gamm4$gam, term = "avgC_capture_day")$data %>%
  dplyr::select(avgC_capture_day = avgC_capture_day, fit = .estimate, se = .se, cil = .lower_ci, ciu = .upper_ci)

dex_cort_plot <- ggplot() +
  geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
  geom_line(data = dex_cort, aes(x = avgC_capture_day, y = fit), color = dex_color, linetype = "solid", size = 0.3) +
  geom_ribbon(data = dex_cort, aes(x = avgC_capture_day, ymin = cil, ymax = ciu), fill = dex_color, alpha = 0.4) +
  labs(x = NULL, y = "Predicted Partial Post-dexamethasone Corticosterone") +
  theme_classic() +
  scale_x_continuous(breaks = seq(0, 30, 10)) +
  xlim(c(13, 30)) +
  theme(panel.grid = element_blank(), axis.text = element_text(size = 12), axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15),
)

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

dex_cort_plot

Combine all plots together

combined <- base_cort_plot | stress_cort_plot | dex_cort_plot  # one row
final_figure1 <- combined + plot_annotation(
    tag_levels = "A",    
    caption = "Average Daily Temperature (°C)"
  ) &
  theme(
    plot.tag = element_text(
      size = 18,
      family = "Helvetica",
      face = "bold"
    ),
    plot.tag.position = c(0.25, 0.98),  # inside the panel (x,y in npc coords)
    plot.caption = element_text(
      hjust = 0.5,
      size = 14,
      family = "Helvetica"
    )
  )
final_figure1

ggsave("TablesFigures/Figure1.png",
       plot = final_figure1, height = 5, width = 14)

Question 2: Do male and female nestlings differ in their corticosterone regulation

Here, I am just using regular linear mixed models to determine if there are any sex differences.

Baseline corticosterone

bleed1_sex <- lmer(log_bleed1_final_cort ~ sex + bleed1_latency_sec + ( 1 | nest_key), data = sex_df)
summary(bleed1_sex)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_bleed1_final_cort ~ sex + bleed1_latency_sec + (1 | nest_key)
##    Data: sex_df
## 
## REML criterion at convergence: 1757.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4701 -0.5800 -0.0154  0.5371  3.4105 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  nest_key (Intercept) 0.6510   0.8069  
##  Residual             0.7263   0.8522  
## Number of obs: 601, groups:  nest_key, 164
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        9.494e-01  1.274e-01 5.389e+02   7.454 3.63e-13 ***
## sexMale            3.989e-02  7.865e-02 5.153e+02   0.507    0.612    
## bleed1_latency_sec 2.010e-03  1.103e-03 5.433e+02   1.822    0.069 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sexMal
## sexMale     -0.282       
## bld1_ltncy_ -0.765 -0.003

simulationOutput = simulateResiduals(bleed1_sex, plot = F)
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0048, p-value = 0.864
## alternative hypothesis: two.sided

Stress-induced corticosterone

bleed2_sex <- lmer(log_bleed2_final_cort ~ sex + log_bleed1_final_cort + ( 1 | nest_key), data = sex_df)
summary(bleed2_sex)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_bleed2_final_cort ~ sex + log_bleed1_final_cort + (1 | nest_key)
##    Data: sex_df
## 
## REML criterion at convergence: 1459.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1207 -0.5130  0.0789  0.5499  2.2702 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  nest_key (Intercept) 0.5201   0.7212  
##  Residual             0.4863   0.6974  
## Number of obs: 572, groups:  nest_key, 158
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)             2.25749    0.08172 295.17103  27.626  < 2e-16 ***
## sexMale                 0.05062    0.06597 475.13834   0.767    0.443    
## log_bleed1_final_cort   0.25653    0.03526 561.00350   7.276 1.16e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sexMal
## sexMale     -0.359       
## lg_bld1_fn_ -0.464 -0.029

simulationOutput = simulateResiduals(bleed2_sex, plot = F)
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0182, p-value = 0.824
## alternative hypothesis: two.sided

Post-dex corticosterone

bleed3_sex <- lmer(log_bleed3_final_cort ~ sex + log_bleed1_final_cort + ( 1 | nest_key), data = sex_df)
summary(bleed3_sex)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_bleed3_final_cort ~ sex + log_bleed1_final_cort + (1 | nest_key)
##    Data: sex_df
## 
## REML criterion at convergence: 957.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5645 -0.5057 -0.0177  0.5134  5.0510 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  nest_key (Intercept) 0.2921   0.5404  
##  Residual             0.3224   0.5678  
## Number of obs: 454, groups:  nest_key, 121
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)             1.79434    0.07105 221.04399  25.256   <2e-16 ***
## sexMale                -0.09935    0.06002 378.72536  -1.655   0.0987 .  
## log_bleed1_final_cort   0.01622    0.03227 446.74111   0.503   0.6154    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sexMal
## sexMale     -0.394       
## lg_bld1_fn_ -0.444 -0.022

simulationOutput = simulateResiduals(bleed3_sex, plot = F)
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97126, p-value = 0.8
## alternative hypothesis: two.sided

Question 3: Does sex mediate the effect of temperature on the corticosterone response in nestlings?

Here, I will be using GAM agains to fit that smooth term effect of temp on cort, to stay consistent with the first set of models.

Baseline corticosterone

bleed1_int_gamm <- gamm4(
  log_bleed1_final_cort ~ sex + s(avgC_capture_day, by = sex) + (bleed1_latency_sec),
  random = ~(1 | nest_key),
  data = sex_df)
summary(bleed1_int_gamm$gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log_bleed1_final_cort ~ sex + s(avgC_capture_day, by = sex) + 
##     (bleed1_latency_sec)
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.966837   0.119479   8.092 3.35e-15 ***
## sexMale            0.050639   0.077758   0.651   0.5151    
## bleed1_latency_sec 0.002023   0.001085   1.865   0.0627 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df     F p-value    
## s(avgC_capture_day):sexFemale 3.216  3.216 17.27  <2e-16 ***
## s(avgC_capture_day):sexMale   3.406  3.406 16.92  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.202   
## lmer.REML = 1719.1  Scale est. = 0.73292   n = 601

summary(bleed1_int_gamm$mer)

## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 1719.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6269 -0.6366  0.0073  0.5611  3.4329 
## 
## Random effects:
##  Groups   Name                          Variance Std.Dev.
##  nest_key (Intercept)                   0.3902   0.6247  
##  Xr.0     s(avgC_capture_day):sexMale   1.4040   1.1849  
##  Xr       s(avgC_capture_day):sexFemale 1.0624   1.0307  
##  Residual                               0.7329   0.8561  
## Number of obs: 601, groups:  nest_key, 164; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                    Estimate Std. Error t value
## X(Intercept)                       0.966837   0.119426   8.096
## XsexMale                           0.050639   0.077617   0.652
## Xbleed1_latency_sec                0.002023   0.001085   1.865
## Xs(avgC_capture_day):sexFemaleFx1 -0.376943   0.314524  -1.198
## Xs(avgC_capture_day):sexMaleFx1   -0.654523   0.351124  -1.864
## 
## Correlation of Fixed Effects:
##             X(Int) XsexMl Xbl1__ X(C__):F
## XsexMale    -0.298                       
## Xbld1_ltnc_ -0.801 -0.004                
## Xs(C__):FF1 -0.035  0.003  0.037         
## Xs(C__):MF1 -0.013  0.002  0.015  0.165

plot(bleed1_int_gamm$gam, pages = 1) # view smooth term

k.check(bleed1_int_gamm$gam) # check the k to see if the model is "wiggly" enough

##                               k'      edf  k-index p-value
## s(avgC_capture_day):sexFemale  9 3.215974 0.819781       0
## s(avgC_capture_day):sexMale    9 3.405542 0.819781       0

par(mfrow = c(2, 2))
gam.check(bleed1_int_gamm$gam) # check residuals to make sure model is normal

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                 k'  edf k-index p-value    
## s(avgC_capture_day):sexFemale 9.00 3.22    0.82  <2e-16 ***
## s(avgC_capture_day):sexMale   9.00 3.41    0.82  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pred_base <- ggpredict(bleed1_int_gamm, 
                         terms = c("avgC_capture_day [all]", "sex [Female, Male]"))

base_int <- ggplot(pred_base, aes(x = x, y = predicted, color = group, fill = group)) +
  geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
  geom_line(size = 1.0) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.4, color = NA) +
  labs(x = NULL, y = "Predicted Baseline Corticosterone") +
  theme_classic() +
  theme(
    panel.grid = element_blank(),
    axis.text  = element_text(size = 12),
    axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15),
    legend.title = element_blank(),
    legend.position = "top") +
  scale_color_manual(values = c(Female = female_color, Male = male_color))+
  scale_fill_manual(values = c(Female = female_color, Male = male_color))+
  scale_x_continuous(breaks = seq(13, 30, 5))
base_int

Stress-induced corticosterone

bleed2_int_gamm <- gamm4(
  log_bleed2_final_cort ~ sex + s(avgC_capture_day, by = sex) + (log_bleed1_final_cort),
  random = ~(1 | nest_key),
  data = sex_df)
summary(bleed2_int_gamm$gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log_bleed2_final_cort ~ sex + s(avgC_capture_day, by = sex) + 
##     (log_bleed1_final_cort)
## 
## Parametric coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.28829    0.07935  28.836  < 2e-16 ***
## sexMale                0.05430    0.06558   0.828    0.408    
## log_bleed1_final_cort  0.23400    0.03558   6.576  1.1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                 edf Ref.df     F p-value    
## s(avgC_capture_day):sexFemale 3.642  3.642 6.012 0.00044 ***
## s(avgC_capture_day):sexMale   2.279  2.279 4.951 0.00707 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.212   
## lmer.REML = 1453.8  Scale est. = 0.48081   n = 572

summary(bleed2_int_gamm$mer)

## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 1453.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1192 -0.4925  0.0871  0.5539  2.4230 
## 
## Random effects:
##  Groups   Name                          Variance Std.Dev.
##  nest_key (Intercept)                   0.4614   0.6792  
##  Xr.0     s(avgC_capture_day):sexMale   0.2509   0.5009  
##  Xr       s(avgC_capture_day):sexFemale 1.4752   1.2146  
##  Residual                               0.4808   0.6934  
## Number of obs: 572, groups:  nest_key, 158; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## X(Intercept)                       2.28829    0.07985  28.657
## XsexMale                           0.05430    0.06543   0.830
## Xlog_bleed1_final_cort             0.23400    0.03589   6.519
## Xs(avgC_capture_day):sexFemaleFx1 -0.48520    0.33575  -1.445
## Xs(avgC_capture_day):sexMaleFx1   -0.24713    0.19585  -1.262
## 
## Correlation of Fixed Effects:
##             X(Int) XsexMl Xl_1__ X(C__):F
## XsexMale    -0.361                       
## Xlg_bld1_f_ -0.490 -0.033                
## Xs(C__):FF1 -0.024  0.002  0.036         
## Xs(C__):MF1 -0.053 -0.007  0.093  0.187

plot(bleed2_int_gamm$gam, pages = 1) # view smooth term

k.check(bleed2_int_gamm$gam) # check the k to see if the model is "wiggly" enough

##                               k'      edf   k-index p-value
## s(avgC_capture_day):sexFemale  9 3.642472 0.6378559       0
## s(avgC_capture_day):sexMale    9 2.279260 0.6378559       0

par(mfrow = c(2, 2))
gam.check(bleed2_int_gamm$gam) # check residuals to make sure model is normal

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                 k'  edf k-index p-value    
## s(avgC_capture_day):sexFemale 9.00 3.64    0.64  <2e-16 ***
## s(avgC_capture_day):sexMale   9.00 2.28    0.64  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pred_stress <- ggpredict(bleed2_int_gamm, 
                       terms = c("avgC_capture_day [all]", "sex [Female, Male]"))



stress_int <- ggplot(pred_stress, aes(x = x, y = predicted, color = group, fill = group)) +
  geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
  geom_line(size = 1.0) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.4, color = NA) +
  labs(x = NULL, y = "Predicted Stress-Induced Corticosterone") +
  theme_classic() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 12),
    axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15),
        legend.title = element_blank(),
        legend.position = "top") +
  scale_color_manual(values = c(Female = female_color, Male = male_color))+
  scale_fill_manual(values = c(Female = female_color, Male = male_color))+
  scale_x_continuous(breaks = seq(13, 30, 5))

stress_int

Post-dex corticosterone

bleed3_int_gamm <- gamm4(
  log_bleed3_final_cort ~ sex + s(avgC_capture_day, by = sex) + (log_bleed1_final_cort),
  random = ~(1 | nest_key),
  data = sex_df)
summary(bleed3_int_gamm$gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log_bleed3_final_cort ~ sex + s(avgC_capture_day, by = sex) + 
##     (log_bleed1_final_cort)
## 
## Parametric coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.812871   0.071565  25.332   <2e-16 ***
## sexMale               -0.094864   0.060123  -1.578    0.115    
## log_bleed1_final_cort -0.001366   0.033306  -0.041    0.967    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf Ref.df     F p-value  
## s(avgC_capture_day):sexFemale 2.15   2.15 2.349  0.1174  
## s(avgC_capture_day):sexMale   2.29   2.29 2.572  0.0979 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0317   
## lmer.REML = 963.09  Scale est. = 0.32267   n = 454

summary(bleed3_int_gamm$mer)

## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 963.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4835 -0.5153 -0.0258  0.5166  5.0105 
## 
## Random effects:
##  Groups   Name                          Variance Std.Dev.
##  nest_key (Intercept)                   0.2779   0.5272  
##  Xr.0     s(avgC_capture_day):sexMale   0.1564   0.3955  
##  Xr       s(avgC_capture_day):sexFemale 0.1199   0.3463  
##  Residual                               0.3227   0.5680  
## Number of obs: 454, groups:  nest_key, 121; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                    Estimate Std. Error t value
## X(Intercept)                       1.812871   0.070667  25.654
## XsexMale                          -0.094864   0.060046  -1.580
## Xlog_bleed1_final_cort            -0.001366   0.033224  -0.041
## Xs(avgC_capture_day):sexFemaleFx1  0.121431   0.146830   0.827
## Xs(avgC_capture_day):sexMaleFx1    0.128565   0.165693   0.776
## 
## Correlation of Fixed Effects:
##             X(Int) XsexMl Xl_1__ X(C__):F
## XsexMale    -0.391                       
## Xlg_bld1_f_ -0.460 -0.031                
## Xs(C__):FF1  0.030  0.003 -0.086         
## Xs(C__):MF1  0.043  0.028 -0.096  0.174

plot(bleed3_int_gamm$gam, pages = 1) # view smooth term

k.check(bleed3_int_gamm$gam) # check the k to see if the model is "wiggly" enough

##                               k'      edf   k-index p-value
## s(avgC_capture_day):sexFemale  9 2.150400 0.6986927       0
## s(avgC_capture_day):sexMale    9 2.290195 0.6986927       0

par(mfrow = c(2, 2))
gam.check(bleed3_int_gamm$gam) # check residuals to make sure model is normal

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                                 k'  edf k-index p-value    
## s(avgC_capture_day):sexFemale 9.00 2.15     0.7  <2e-16 ***
## s(avgC_capture_day):sexMale   9.00 2.29     0.7  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pred_dex <- ggpredict(bleed3_int_gamm, 
                         terms = c("avgC_capture_day [all]", "sex [Female, Male]"))


dex_int <- ggplot(pred_dex, aes(x = x, y = predicted, color = group, fill = group)) +
  geom_abline(slope = 0, intercept = 0, linetype = "61", color = "gray75", size = 0.3) +
  geom_line(size = 1.0) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.4, color = NA) +
  labs(x = NULL, y = "Predicted Post-dexamethasone Corticosterone") +
  theme_classic() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 12),
    axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15, 15, 15, 15),
        legend.title = element_blank(),
        legend.position = "top") +
  scale_color_manual(values = c(Female = female_color, Male = male_color))+
  scale_fill_manual(values = c(Female = female_color, Male = male_color))+
  scale_x_continuous(breaks = seq(13, 30, 5))
dex_int

Combine all together

combined2 <- base_int | stress_int | dex_int  # one row

final_figure2 <- combined2 + plot_annotation(
    tag_levels = "A",    
    caption = "Average Daily Temperature (°C)"
  ) &
  theme(
    plot.tag = element_text(
      size = 18,
      family = "Helvetica",
      face = "bold"
    ),
    plot.tag.position = c(0.25, 0.85),  
    plot.caption = element_text(
      hjust = 0.5,
      size = 14,
      family = "Helvetica"
    ))
  
  final_figure2

ggsave("TablesFigures/Figure2.png",
       plot = final_figure2, height = 5, width = 14)

Female / Male Comparision Plots

Now, creating a figure to compare the smooths in the interaction models, based on this tutorial: this tutorial:https://fromthebottomoftheheap.net/2017/10/11/difference-splines-i/

Run function from above site:

smooth_diff <- function(model, newdata, f1, f2, var, alpha = 0.05,
                        unconditional = FALSE) {
  fac <- as.character(newdata[[var]]) # added this to convert sex to a character to allow for comparision
  xp <- predict(model, newdata = newdata, type = 'lpmatrix')
  c1 <- grepl(f1, colnames(xp))
  c2 <- grepl(f2, colnames(xp))
  r1 <- fac == f1 # modified this
  r2 <- fac == f2
  ## difference rows of xp for data from comparison
  X <- xp[r1, ] - xp[r2, ]
  ## zero out cols of X related to splines for other lochs
  X[, ! (c1 | c2)] <- 0
  ## zero out the parametric cols
  X[, !grepl('^s\\(', colnames(xp))] <- 0
  dif <- X %*% coef(model)
  se <- sqrt(rowSums((X %*% vcov(model, unconditional = unconditional)) * X))
  crit <- qt(alpha/2, df.residual(model), lower.tail = FALSE)
  upr <- dif + (crit * se)
  lwr <- dif - (crit * se)
  data.frame(pair = paste(f1, f2, sep = '-'),
             diff = dif,
             se = se,
             upper = upr,
             lower = lwr)
}

Create dataframe to be fed into the function:

pdat_base <- expand.grid(
  avgC_capture_day = seq(13.5, 28.5, length = 400),
  sex = c("Female", "Male"),
  bleed1_latency_sec = mean(sex_df$bleed1_latency_sec, na.rm = TRUE)
)

pdat_stress_dex <- expand.grid(
  avgC_capture_day = seq(13.5, 28.5, length = 400),
  sex = c("Female", "Male"),
  log_bleed1_final_cort = mean(sex_df$log_bleed1_final_cort, na.rm = TRUE)
)

comp_base <- smooth_diff(bleed1_int_gamm$gam, pdat_base, 'Female', 'Male', 'sex')
comp_base$avgC_capture_day <- pdat_base$avgC_capture_day[pdat_base$sex == "Female"]
comp_stress <- smooth_diff(bleed2_int_gamm$gam, pdat_stress_dex, 'Female', 'Male', 'sex')
comp_stress$avgC_capture_day <- pdat_stress_dex$avgC_capture_day[pdat_stress_dex$sex == "Female"]
comp_dex <- smooth_diff(bleed3_int_gamm$gam, pdat_stress_dex, 'Female', 'Male', 'sex')
comp_dex$avgC_capture_day <- pdat_stress_dex$avgC_capture_day[pdat_stress_dex$sex == "Female"]

Create plots:

Baseline Corticosterone

base_diff <- ggplot(comp_base, aes(x = avgC_capture_day, y = diff)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = purple_color) +
  geom_line(size = 1, color = purple_color) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_cartesian(ylim = c(-1, 1)) +
  labs(x = NULL, y = 'Estimated Difference in Baseline Corticosterone \n(Female - Male)') +
  theme_classic() +
  theme(panel.grid = element_blank(), axis.text = element_text(size = 12),axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15,15,15,15))
base_diff

## Stress-induced corticosterone

stress_diff <- ggplot(comp_stress, aes(x = avgC_capture_day, y = diff)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.4, fill = purple_color) +
  geom_line(size = 1, color = purple_color) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_cartesian(ylim = c(-1, 1)) +
  labs(x = NULL, y = 'Difference in Stress-Induced Corticosterone\n(Female - Male)') +
  theme_classic() +
  theme(panel.grid = element_blank(), axis.text = element_text(size = 12),axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15,15,15,15))
  
stress_diff

Post-dex corticosterone

dex_diff <- ggplot(comp_dex, aes(x = avgC_capture_day, y = diff)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.4, fill = purple_color) +
  geom_line(size = 1, color = purple_color) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  coord_cartesian(ylim = c(-1, 1)) +
  labs(x = NULL, y = 'Difference in Post-dexamethasone Corticosterone\n(Female - Male)') +
  theme_classic() +
  theme(panel.grid = element_blank(), axis.text = element_text(size = 12),axis.title = element_text(size = 14, family = "Helvetica"), plot.margin = margin(15,15,15,15))
dex_diff

Create one figure with all combined:

combined3 <- base_diff | stress_diff | dex_diff  # one row
final_figure3 <- combined3 +
  plot_annotation(
    tag_levels = "A",    
    caption = "Average Daily Temperature (°C)"
  ) &
  theme(
    plot.tag = element_text(
      size = 18,
      family = "Helvetica",
      face = "bold"
    ),
    plot.tag.position = c(0.25, 0.98), 
    plot.caption = element_text(
      hjust = 0.5,
      size = 14,
      family = "Helvetica"
    )
  )

final_figure3

ggsave("TablesFigures/Figure3.png",
       plot = final_figure3, height = 5, width = 14)

nestling-temp-cort-final

Erica Griffin

2025-11-23

Environmental effects on corticosterone in nestling tree swallows

Last updated: November 17, 2025

Abstract:

Research Questions:

Hypotheses:

Predictions:

Question 1: How does temperature influence nestling corticosterone levels?

Baseline corticosterone

Stress-induced corticosterone

Dex corticosterone

Question 2: Do male and female nestlings differ in their corticosterone regulation

Baseline corticosterone

Stress-induced corticosterone

Post-dex corticosterone

Question 3: Does sex mediate the effect of temperature on the corticosterone response in nestlings?

Baseline corticosterone

Stress-induced corticosterone

Post-dex corticosterone

Female / Male Comparision Plots

Baseline Corticosterone

Post-dex corticosterone